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1 Stable and Asymptotocally Stable Com- 
pact Schemes 

Recently, higher-order compact schemes have seen increasing use in the DNS 
(Direct Numerical Simulations) of the Navier-Stokes equations. Although 
they do not have the spatial resolution of Spectral methods, they offer signif- 
icant increases in accuracy over conventional second order methods. They can 
be used on any smooth grid, and do not have an overly restrictive CFL depen- 
dence as compared with the 0(N ~ 2 ) CFL dependence observed in Chebyshev 
Spectral methods on finite domains. In addition, they are generally more ro- 
bust and less costly than Spectral methods. The issue of the relative cost 
of higher-order schemes (accuracy weighted against physical and numerical 
cost) is a far more complex issue, depending ultimately on what features of 
the solution are sought and how accurately they must be resolved. In any 
event, the further development of the underlying stability theory of these 
schemes is important. 

It turns out that this schemes are very sensitive to boundary treatments. 
In particular all of the boundary conditions , currently used , allow non 
physical time growth of the solution. Recently, the stability characteristics 
of various compact fourth- and sixth-order spatial operators were assessed in 
reference [1], using the theory of Gustafsson, Kreiss and Sundstrom (G-K-S) 
for the semi-discrete Initial Boundary Value Problem (IBVP). The results 
were then generalized to the fully discrete case with Runge-Kutta time ad- 
vancement using a recently developed theory by Kreiss. In all cases, favorable 
comparisons were obtained between G-K-S theory, eigenvalue determination, 
and numerical simulation. The conventional definition of stability is then 
sharpened to include only those spatial discretizations that are asymptoti- 
cally stable (bounded, Left Half-Plane eigenvalues). It is shown that many 
of the higher-order schemes which are G-K-S stable are not asymptotically 
stable. It was concluded that in practical calculations, only those schemes 
which satisfied both definitions of stability were of any great usefulness. 

It was shown in the above work of Carpenter et al. that conventional (op- 
timal) finite difference closures at the boundaries of order greater than four 
are not G-K-S (or asymptotically) stable. Since fifth-order boundary clo- 
sures possessing both stability properties were needed for sixth-order inner 
schemes, an alternate method for closing the boundaries was sought. The so- 
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lution was to parametrize the fifth-order difference formula at several points 
at each end of the spatial domain, thereby creating adjustable coefficients in 
the spatial operator. The asymptotic properties of the operator were estab- 
lished by the numerical determination of the eigenvalue spectrum, and the 
parameters were then adjusted until the desired spectrum was obtained. The 
resulting scheme was then tested for G-K-S stability, and if stable, satisfied 
both the desired criteria for a numerical discretization. 

Several technical difficulties were encountered in trying to determine sta- 
ble formulations in this manner. In general, a large number of free parameters 
were needed to find a combination which resulted in a stable formulation. 
This results from trying to achieve a high-order discretizations at the in- 
flow boundary where the stencils are dramatically downwind, and mostly 
unstable. Although a stable closure condition was found for the sixth-order 
compact scheme, (5 2 ,5 2 — 6 — 5 2 ,5 2 ) it was apparent that if schemes of higher 
accuracy were to be obtained, a systematic procedure was required to con- 
strain the parameter space over which the search was performed. Another 
difficulty was that the numerical eigenvalue determination did not yield the 
exact eigenvalues of the spatial operator, but rather depended on numerical 
round off and the condition number of the resulting spatial operator. This 
was not found to be a significant problem for the schemes determined in the 
study, but it was found that many of the high-order schemes were not well 
conditioned. 

The fundamental difficulty with determining a spatial operator based 
on the results from an eigenvalue analysis, is that it uses as a basis for the 
method the the spatial matrix resulting from discretization of the scalar wave 
equation U t + aU x — 0. While G-K-S stability of a discretization on 
the scalar wave equation implies G-K-S stability on a system of 
hyperbolic equations, (if the boundary conditions are imposed in 
characteristic form) [2], the same is not in general true for asymp- 
totic stability. Therefore, there is no guarantee that the numerical scheme 
determined in this manner will be stable for an arbitrary hyperbolic problem. 
. An obvious remedy for the analysis presented by Carpenter et al.[l] would 
to have used the system not the scalar eigenvalue determination as a basis 
for devising stable closure formula. This would further constrain an already 
difficult search procedure to isolate the parameters at the boundaries which 
would produce an strictly asymptotically stable scheme. An altogether dif- 
ferent procedure must be used if an arbitrarily high-order scheme is sought. 
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The approach of devising suitable boundary closures and then testing 
them with various stability techniques (such as finding the norm) is entirely 
the wrong approach when dealing with high-order methods. Very seldom 
are high-order boundary closures stable, making them difficult to isolate. 
An alternative approach is to begin with a norm which satisfies all the sta- 
bility criteria for the hyperbolic system, and look for the boundary closure 
forms which will match the norm exactly. This method was used recently 
by Strand[4] to isolate stable boundary closure schemes for the explicit cen- 
tral fourth- and sixth-order schemes. The norm used was an energy norm 
mimicking the norm for the differential equations. Further research should 
be devoted to BC for high order schemes in order to make sure that the 
results obtained are reliable. The compact 4-th order and sixth order finite 
difference scheme had been incorporated into the a code to simulate flow past 
circular cylinders. This code will serve as a verification of the full spectral 
codes. A detailed stability analysis by Carprnter (from the fluid Mechan- 
ics Division) and Gottlieb gave analytic conditions for stability as well as 
asymptotic stability.This had been incorporated in the code in form of stable 
boundary conditions. 

Effects of the cylinder rotations has been studied. The results differ from 
the known theoretical results. We are in the middle of analyzing the results. 

A detailed analysis of the effects of the heating of the cylinder on the 
shedding frequency had been studied using the above schemes. It has been 
found that the shedding frequency decreases when the wire was heated. Ex- 
perimental work is being carried out to affirm this result. This is carried out 
by Eric Voth in conjuction with D. Rudy from the Fluid Mechanics Division. 
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2 Wavelets 

A major effort to adapt wavelets to the solution of PDE’s is under investiga- 
tion. It has been found by L. Jameson ( agraduate student in the program) 
that using wavelets as a basis function for differentiation is equivalent to the 
use of finite difference schemes. The result is suppose to give a clue of how to 
implement boundary conditions. We attach a paper by him on the subject. 
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0,1 Introduction 


The numerical solution of a partial differential equation requires an easily 
manipulated spatial approximation to the derivative of the unknown func- 
tion as well as some method to march forward in time. In general one starts 
from given values of the unknown function, then a finite dimensional ap- 
proximation, based on those values, is constructed. This approximation is 
differentiated and the result are read at the gridpoints. For example, in the 
psuedospectral Chebyshev method for the disretization of the equation 

dU(x,t) _ dF[U(x,t)] 
dt dx 

One assumes that at a given time the values of U(xj,t) are given for some 
points xj — cos ( -jy ) , (j = 0, TV). Then one constructs the interpolation 
polynomial through those points and differentiate this polynomial to get ap- 
proximate values for at the point Xj. this procedure can be viewed 

as a transformation from N given values (of the function) to new N values 
(approximating the derivative. This is the Chebyshev Differentiation Matrix . 
The numerical algorithm therefore is simple and the boundary conditions can 
be easily applied. It is natural to ask whether one gain by using wavelets 
instead of Chebyshev polynomials. . Since wavelets are well localized func- 
tions it is reasonable to conjecture that they might represent steep gradients 
or the development of a shock with a relatively small number of terms, con- 
sider a periodic function f(x) given on an equally spaced mesh. Expand it 
in wavelet expansion and use the derivative of this expansion as an approxi- 
mation to the derivative of /(x). Lee Jameson, a student of D. Gottlieb has 
recently proved that approximation of a periodic function f(x) in a wavelet 
basis and the differention of this approximation yields nothing more than a 
finite difference approximation to a derivative. The following is an outline of 
the proof: 

i) Given a periodic function f(x) let s represent the periodic scaling func- 
tion coefficients of this function on the finest scale. This requires approx- 
imating the inner product of f(x) with the scaling function on the finest 
scale. The matrix representation of this approximation is circulent in form: 
C : / — ► s , where / represents f(x) sampled on an evenly spaced grid. 

ii) Let D be the mapping from the scaling function coefficients of f(x) to 
the set of scaling function coefficients that represents the derivative of f(x ): 
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D : s — ► s. Since f(x) is periodic then the matrix form of D is circulent in 
form. 

iii) All circulent martrices of the same size commute, therefore we can 
apply the operator D directly to /. The operator D has the effect of a finite 
difference operator, and the proof will be complete. 

Therefore, the effect of first approximating in a wavelet basis, then differ- 
entiating in this basis and finally converting back from the wavelet basis to 
the original function is equal to applying the appropriate finite difference op- 
erator directly to the equally spaced sampled values of the original function 
f(x). The proof provides an insight into the possibility of using wavelets for 
solutions of PDE’s. 

Wavelets while not more than known finite difference schemes 
can provide a mechnism for automatic adaptation of the mesh. 

Since the proof is unpublished we will bring it here in some detail. This 
proof contains five sections: i) The first is the introduction outlining the pre- 
sentation. ii) The second introduces scaling functions and wavelets, iii) The 
third discusses the approximation of a periodic function by scaling functions 
on the finest scale, iv) The fourth is concerned with the derivative of the scal- 
ing function and wavelet approximation of a function, v) The fifth concludes 
with a statement of the thesis that spatial wavelet approximations provide 
nothing more than finite difference methods do for the numerical solution of 
partial differential equations. 

0.2 Definition of Wavelets 

Wavelets have been precisely defined in many places [Daubechie], [Strang], 
[Beylkin], and others. The following outlines the most prominent properties 
of wavelets. 

Begin with two sets of coefficients of length T, [Daubechie] H = { A: } * = o , 
G = {gk}kZo called quadrature mirror filters, i.e., H and G are related by 
Qk = (— 1 for k = 0, ...,L — 1, which completely determine, along with 

the additional restriction of normalization, 

<f)(x)dx = 1, 

the mother scaling function <f>(x) and the mother wavelet ^(x), respectively, 
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by the following relations: 


L - 1 

<f >( x ) = ^2 hk<f>(2x - k), 

Jt=0 


and 


L—l 


H x ) = ]C gk<f>( 2x ~ k )- 

k = 0 


For the remainder of this paper the following notation will adopted: <f>i(x) 
and ipi(x) will denote the mother scaling function and mother wavelet, re- 
spectively, at scale j and location k, i.e., 


4M = - k), 


and 

~ — k ). 

A few of the ramifications of the above definitions are, first of all, that 
the wavelet ift(x) has M vanishing moments 

/ oo 

ift(x)x m dx — 0 

■CO 

for m = 0 , M — 1, where the number of coefficients in H and G is equal to 
twice the number of vanishing moments, L = 2 M (this is true for the usual 
Daubechie wavelets only). Second, define Vj and Wj as linear span of (f){ and 
if> J k over all location parameters k with the scale j fixed, i.e., 

Vj = span k <f> J h (x) y 


and 

Wj = span k r() J k (x). 

These definitions lead to, 

... C^CVoC C 

n^ = w, 

jez 
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U v i = i 2 («). 

)6Z 

^+i = Vi© Wj. 

and finally, 

i ! (R) = ® HO. 

j€Z 

These are the essential definitions and properties of scaling functions and 
wavelets. For more discussion and details see [Daubechie], [Mallat], and 
[Strang]. 

0.3 Approximating in Wavelet Bases 

Scaling functions and wavelets were introduced in the previous section. As 
noted, Vj is the space spanned by <j>k{x) over all k . Without loss of generality, 
let scale j = 0 be the finest scale. Then, for example, Vo = V\ © W\. The 
approximation of an arbitrary periodic function f(x) begins by projecting 
f(x) onto each basis function 4>^(x) at the finest scale: 

/ OQ 

f(x)<f>° k (x)dx. 

-OO 

The approximation properties of scaling functions are determined by the 
number of vanishing moments of the associated wavelet: if the mother wavelet 
has M vanishing moments then the polynomials 1, x, x Af_1 are lin- 
ear combinations of the translates of the mother scaling function <f>(x — k ) 
[Strang]. Furthermore, smooth functions can be approximated with error 
C )(h M ) [Strang], where h represents the grid size. 

Once the function /(x ) has been approximated on the finest scale, j = 0, 
then the coefficients s® can be decomposed into coefficients at scales that are 
twice as course at each decomposition using the following equations [Mallat]: 

n—2M 

= X] ^" 5 n+2A:-2 

n=l 

n=2M 

= 9n S n+2k-2‘> 

n= 1 
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where M is the number of vanishing moments of the wavelet and {h} and 
{</} are the quadrature mirror filters defined in the previous section. 

To restate, first the function f(x) is approximated at the finest scale j = 0 
with an error of order M to the coefficients s® then the coefficients at more 
coarse scales are found by the above pyramid-like decomposition. 

A second method suggested by Beylkin et. al. [Beylkin] is to approximate 
the integral of each scaling function coefficient and each wavelet coefficient 
directly from the integral, i.e., by an appropriate quadrature formula approx- 
imating the following integrals: 

/ oo 

f(x)<f>{(x)dx, 

-OO 

/ OO 

f(x)tl>l(x)dx. 

-oo 

In this paper the first method will be used so that all approximations 
will be made at the finest scale. Furthermore, all work will be done with the 
usual Daubechie wavelets. For wavelets supported on [0, 3M] see appendix 
1. 


0.4 Quadrature Formula for Scaling Function 

The scaling function coefficients of a function f(x) on the finest scale are 
calculated exactly by, 

/ OO 

f(x)4>{x — k)dx. 

-oo 

For a numerical calculation, however, one must work with an estimate of the 
above coefficients, sjj!, i.e., a suitable quadrature formula is needed. 

Recall from the previously stated approximation properties of scaling 
functions that if the associated wavelet has M vanishing moments then one 
can represent polynomials up to order M — 1 exactly by translations of the 
scaling function <f>(x ). Therefore, for f(x) equal to a polynomial up to order 
M — 1 the scaling function coefficients, s£, can be found exactly. Conse- 
quently, there exist a set of coefficients {c/}^ 1 such that 

M—l 

f(y + k)(f>(y)dy = ^ c,f(l + k), 

l=o 
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where f(x) can be a polynomial up to degree M — 1, and the above integral 
is a shifted version of f(x)<f)(x — k)dx. More simply, the coefficients 

o 1 can found [Beylkin] by solving the following linear system: 


M—\ 


E '“<* = 


/ oo 

-OO 


for m = 0, 1, Af — 1. 

Note that the above quadrature formula will yield an estimate of s° k with 
error of order M. Also, note that the derivation of the coefficients of the 
quadrature formula depend only on the moments of the scaling function 
f^° CG x m cf>(x)dx, 


0.5 Example with D$ 

The ideas in this paper are quite simple and are probably best illustrated by 
an example. The example will be for the Daubechie [Daubechie] wavelet D$. 
The objective is to derive the matrix form of the mapping from evenly-spaced 
samples of a periodic function f(x) to the scaling function coefficients on the 
finest scale Comparable results for the wavelets Z) 4 and Dg are presented 
in appendix 2. 

Recall, first of all, that in the previous subsection the coefficients {q}^ 1 
were determined from the moments of the scaling function. Therefore, the 
scaling function moments must first be calculated. 

Let Mi be the l — ih moment of the scaling function <f>(x): 

Mi = J <f>(x)x t dx , 

and let /q be the / — th moment of the filter h^: 

m = ^2 kl hk- 

k 

Recall that it is required that <f>{x ) be normalized: 

M° ~ J = h 

Also, by integrating the definition of <f>(x) the following results: 

J (f>(x)dx = Y,hkj <f>(2x — k)dx . 

k 
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Let y = 2x — k to get, 


which implies, 


= \ ll h k J <f>iv) d y . 


/*o = — 2. 

it 


The /x/ for / > 0 can be found by straight-forward calculation. The M\ for 
l > 0 can be found from [see appendix 2] 


E ( 7 ) 


For the examples presented here moments up through the third moment are 
needed: M x = \ fi u M 2 = |((//i) 2 + ^ 2 ), and M 3 = ^((/ii) 3 + + 2 p 3 ). 

After the appropriate moments are found, the coefficients {c/}£fo 1 can t> e 
found from 

m- 1 f 

£ l m ci = / x m <j>{x)dx 

1=0 d 

for m = 0, 1 , M — 1 . Specifically, for the De wavelet the linear system in 
matrix form appears as, 


1 1 1 
0 1 2 
0 1 4 


ci = M x 
C2 {M2 


which has the solution cq = .1080, Ci = .9667, and c 2 = —.0746. In tabular 

i Mi m Cj 

form, the complete results for Dr are, , JL- Recall 

1 .8174 1.6348 .9667 

2 .6681 1.3363 -.0746 

that the quadrature formula has the form, 

M— 1 1 

4= £«/('+ *) + 0(t 7) m , 

/=0 iV 


where N is the number of points in the grid. If the function f(x) is periodic 
then in matrix notation the above operation is s = Cf where C for D e and 
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a grid of 6 points appears as, 


/ .108 .967 -.075 0 0 0 \ 

0 .108 .967 -.075 0 0 

0 0 .108 .967 -.075 0 

0 0 0 .108 .967 -.075 

-.075 0 0 0 .108 .967 

V .967 -.075 0 0 0 .108 

The important point here is that the above matrix is circulent [Strange’s 
book] in form. This is the most important observation in this paper, be- 
cause all circulent matrices can be diagonalized by the Fourier matrix, i.e., 
all circulent matrices of the same dimensions have the same eigenvectors 
and therefore they commute. The importance of this property will become 
apparent after the wavelet derivative is discussed in the next section. 

0.6 Derivative based on Wavelets 

In the previous section the mapping from evenly-spaced samples of a periodic 
function, /(x), to the scaling function coefficients on the finest scale, sj!, was 
discussed. The mapping is nothing more than a quadrature formula which 
is exact for /(x) equal to a polynomial up to order M — 1, where M is the 
number of vanishing moments of the wavelet. The question now is what is 
the mapping from to the coefficients of the derivative of /(x), i.e., the 
scaling function coeffiecients, i°, of f'(x). The answer is provided by Beylkin 
[Beylkin], and is presented in the following subsection. This section of the 
paper is organized as follows: i) Beylkin’s results on derivative projections 
will be presented, ii) It will be argued that one need only consider the 
derivative mapping acting on the scaling function coefficients at the finest 
scale, iii) The similarity between the coefficients derived by Beylkin to finite 
difference approximations to the derivative will be presented. 

0.7 Wavelet Coefficients of the Derivative 

An arbitrary wavelet expansion of a function might contain wavelet coeffi- 
cients and scaling coefficients at many scales. Beylkin derives the projection 
coefficients that map from scaling function coefficients and wavelet function 
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coefficients at a given scale to the derivative scaling function coefficients and 
wavelet function coefficients at the same scale. The matrix elements of these 
projections are computed from, 

/ oo . , 

tp(2~ 3 x — i)ip(2~ 3 x — l)dx, 

-oo 

/ °° . , 

V>(2- J x - i)<f>(2~ 3 x - l)dx, 

-OO 

/ oo . , 

<j>( 2~ 3 x — i)ip(2~ 3 x — l)dx , 

-oo 

/ oo , 

4>(2~ 3 x - i)<f>( 2~ 3 x - l)dx. 

-oo 

It is important to note that these projections are all at the same scale j, and 
that projections across different scales appear to be too complicated to yield 
closed-form solutions. The above projections, however, yield equations that 
are simple to work with [Beylkin]. 

It will be argued in the next section that in order to understand the 
numerical properties of the above projections it is only necessary to consider 

/ oo A 

(f)(x - /)— <f>(x)dx, 

-oo uX 

for / E Z. 


0.8 Derivative of Scaling Function Only 

Before beginning the main argument of this subsection some new notation 
will be introduced. The vectors h and g contain the coefficients of the 
quadrature mirror filters which define the mother scaling function and mother 
wavelet, respectively. Define the unitary projection matrix P as, 


PnxN = 


' h 0 

0 h 

0 0 
9 0 
0 9 

0 0 


0 

0 

0 

0 


5 
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where the matrix subscripts denote the size of the matrix. Of course, P is 
nothing more than a matrix with the vectors h and g placed in its rows with 
the vector shifted two places to the right with each subsequent row. Also, 
let P contain the scaling function coefficients at scale j. P is, therefore, the 
matrix that maps P onto s 5-1 " 1 and d 3+1 . Note that the vectors at scale j + 1 
are half as long as the vectors as scale j: 



Let the four matrices R , A, 5 , C be the derivative projections defined 
in the previous subsection, i.e., Beylkin’s coefficients, are R *-* r{j , A <-» a,j, 
B f 3 t j , and C «-* 7 ,j. Explicitly, the mappings are 


R : P 


5 , 


-k -J 

A : d } d , 

B : P 2 , 

C : d 3 s . 

That is, if P and d? denote the scaling and wavelet coefficients of a function 

-*j 

at scale j then s and d denote the scaling and wavelet coefficients of the 
derivative of the function at the same scale. 

For further illustration, suppose that a periodic function has been approx- 
imated on a grid with 16 scaling function coefficients. One application of the 
above defined matrix Piexi6 on the vector followed by the application of 
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the matrix Pg x $ on the vector s 1 would appear as, 



In the above decomposition there are three ways to represent exactly the same 
information: i) All information is at scale 0, i.e., use only the coefficients s ® 
for i = 1, ..., 16. In this case the derivative coefficients would be found by 
applying the above defined matrix [Beylkin] i? 16xl6 . ii) All information is 
at scale 1, i.e., use the coefficients s } and d ] for i = In this case 

all four of the above defined matrices i? 8x8 , A 8x8 , Z? 8x8 , and C 8x8> but the 
application of these four matrices is exactly the same as applying Riexie at 
scale 0 as is scenario (i). iii) This third scneario is the most unwieldy. The 
information is contained in two scales: the eight coefficients at scale 1, s } for 
i = and eight coefficients at scale 2, s? and d? for i = The 

difficulty here is the projection across scales. That is, how does one project 
the derivative of the wavelet at scale 2 onto the scaling function coefficients 
at scale 1. An attempt to calculate this projection has been made by this 
author but without success. One can, of course, approximate this projection 
but this is not very accurate and not elegant. 

Recall that the main argument of this subsection is to illustrate that it is 
only necessary to take the derivative of a wavelet expansion on the finest scale, 
j = 0. First note that regardless of how the the information is represented in 
each of the above three scenarios there are always 16 degrees-of-freedom, i.e., 
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it does not matter which 16 parameters are used to represent the function and 
its derivative. More explicitly, in the first scenario the derivative coefficients 

■hO 

s are calculated by applying /? 16xl6 to s°: 

2° __ o ~o 

s — ^ 16 x 16 ' S . 


In the second scenario, however, one must first apply P 16xl6 to 5° to get the 
scaling and wavelet coefficients at scale 1. The derivative coefficients at scale 
1 are then calculated by applying /) 16xl6 , where 


Avx/v 


Rn/2xN/ 2 Cfsjf 2XN/2 
Bn/2xN/ 2 Atf/2xN/2 


To clarify, in scenario 2 the following operations are performed: 



= D 


16x16 


16X16 




However, scenarios 1 and 2 are exactly the same since 
Rn*N = 1/2 (PffxN * &NX.N ■ Pnxn)- 

In summary, an attempt has been made to illustrate that the derivative 
coefficients of a scaling and wavelet expansion can be calculated at any scale. 
The goal for this author is to understand exactly what wavelets are and what 
they are doing, therefore, scale 0 provides the clearest scenario in which to 
work without sacrificing essential properties of wavelets. 

Given, now, that it is sufficient to work on scale 0 to understand exactly 
what the wavelet derivative does, one must understand the ramifications of 
applying the matrix R^ x n to the vector s°. In the next subsection the 
similarity between the above defined matrix R and finite difference formulas 
for taking the derivative will be explored. 


0.9 Wavelet Derivatives and Finite Difference 

As the previous subsection attempted to illustrate, the essential properties 
of the wavelet derivative are contained in the matrix R. It was suprising, 
at least to this author, that the elements of the matrix R could differentiate 
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not only the vector s but also the equally spaced samples of a function /(x), 
i.e., the matrix R displays finite difference properties. First of all, it is useful 
to simply note the similarity between finite difference coefficients and the 
elements of the matrix R. The following is a table of centered finite difference 
coe fficients and the order of accuracy of the approximation to the derivative: 


Order of Accuracy 

Coefficients 

2 

-5 o i 

4 

2 2 

i-2 0 2_i 

6 

8 

12 3 u 3 12 

13 3 n 3 3 1 

„ 60 20 4 U 4 20 60 

1 41 4 n 4 1 4 1 

280 1Q5-5 5 U 5 5105 280 


Recall that the elements of the matrix R calculated by Beylkin provide the 
transformation from scaling function coefficients of a function to the scaling 
function coefficients of the derivative of the same function. The coefficients 
for the D 2 and D 4 wavelet expansions are exactly the same as the coefficients 
for the 2-nd and 4-th order centered finite difference formulas. The coeffi- 
cients for the Z) 6 and Dg wavelets are not exactly the same but are essentially 
the same (in a finite difference sense). The similarity can be seen clearly by 
plotting the finite difference coefficients and the wavelet coefficients on the 
same plot. The wavelet coefficients are, 


Wavelet 


Convolution Coefficients 

D 2 



Hoi 
2 u 2 

D 4 

A 
D 8 

1 

16 

J 2 q 2 __1_ 

M \ll n 2Z2 2 53 16 1 

365 365 U 365 365 1095 2920 

1664 2645 128 1 

2920 

39296 

1095 

76113 

49553 

396424 

49553 1189272 743295 1189272 


If the above coefficients are treated as finite-difference coefficients then 
it would be nice to know the accuracy. To establish the finite-difference 
accuracy of the coefficients calculated by Beylkin note that a centered-finite- 
difference derivative approximation with 2 K coefficients, can be 

written 

K 

f( x i) ~ J2 a k{fj+k - fj-k)- 

k=l 

If the above equation is exact for f(x) = x r for r = but not for 

r = N + 1 then the equation is said to be N-th order accurate. Therefore, 


13 



one must check to see if 


rx j 1 = “*;-*)> 

k = 1 

when /(x) = x T . To simplify, one can let Xj = j and check the following: 

rj'- , = E«*(C> + *)'- O' - *D- 

fc=l 

Now, treating the coefficients derived by Beylkin as nothing more than finite- 
difference coefficients one can check the accuracy. The following table con- 
tai ns the results of applying the Beylkin coe fficients to various polynomials: 

Wavelet Exact for But not for 


D 2 


p 

D 4 

X 4 

x 5 

D 6 

x 6 

x 7 

D 8 

00 

H 

x 9 

D 10 

x 10 

X 11 


The pattern in the above table is obvious and leads to the following 
conjecture: the cofficients derived by Beylkin which map scaling function 
coefficients of a function to the scaling function coefficients of the derivative 
of the function for the Daubechie wavelet Z? 2 M can differentiate, exactly, a 
polynomial of degree 2M when applied to the samples of the polynomial in 
a finite-difference sense. 

This concludes the second important observation of this paper (the first 
concerned approximating a function by scaling-function coefficients). The 
previously defined matrix R has as its elements the coefficients which display 
this finite-difference quality, and if the original function f(x) is periodic then 
the matrix R is circulent in form. The following concluding section of this 
paper should unify the presentation. 

0.10 Conclusion 

The two important sections of this paper are sections 3.4.3 and 3.4.44. In 
section 3.4.3 it was established that if given a periodic function f(x) then 
the scaling function coefficients s of the function at the finest scale can be 
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approximated by a quadrature formula which in matrix form, 




yields a circulent matrix C. In section 4 it was noted that the coefficients 
derived by Beylkin which map the scaling function coefficients of a periodic 
fuction to the scaling function coefficients of the derivative of the function is 
also circulent in form when written in matrix notation, 

s' — Ds \ 

Furthermore, the matrix D can differentiate (in a finite-difference sense) 
polynomials exactly up to the order of the wavelet. Now, combine the results 
of sections 3 and 4 to get the following relation: 

/' = C- 1 DC f. 

Throughout the paper it has been noted that C and D are circulent in form 
when f(x) is periodic. Circulent matrices of the same dimensions can, how- 
ever, be diagonalized by the same matrix, the Fourier matrix of appropriate 
dimensions, and this implies that all circulent matrices of the same dimen- 
sions commute. Therefore, the previous relation simply becomes, 

/' = Df, 

but when D is applied to the samples of a function it acts as a finite-difference 
operator. In conclusion, when wavelets are used to solve partial differen- 
tial equations numerically they appear to provide nothing more than finite- 
difference methods provide. 
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0.12 Appendix 1 

0.13 Wavelets Supported on (0,3M) 

In this appendix our wavelets are supported on [0,3A/] where M is the num- 
ber of vanishing moments of the wavelet. These are not the usual Daubechie 
wavelets, but for these wavelets the scaling function coefficients of a periodic 
function f(x) can be approximated with error of order M simply by sampling 
f(x) at the correct location. 

To begin, assume that there exist a unique tm , fixed for a fixed number 
of vanishing moments, M , of the wavelet, such that 

J <j>(x + T\f)x m dx = 0 

for m = 1,2, ...,M — 1. Furthermore, recall the definition of the scaling 
function coefficient and expand the integrand fix) in a Taylor series about 

*o tffi = /’(*„)): 

s° k = j f(x)<f>(x — k)dx = 

fo J 4>(x — k)dx + f' Q J (x — x 0 )(j){x — k)dx + /„ J (x — x 0 ) 2 <f>(x — k)dx + .... 

Now, shift the variable of integration by y = x — r — k, and choose the point 
of expansion, a: 0 , to be r + k to get, 

s° k = 

fi T + k ) J 4>(y + r)dy + f'(y + k) J y4>(y + r)dy + f"(y + k) J i / 2 <f>(y + T)dy + .... 
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Now, rename r as tm and the above integrals are of the form. 


J <f>(x + r M )x m dx = 0 , 

and therefore vanish for rn = 1 — 1 leading to, 

s °k = f{r M + k) + / (M) (tm + k) J y M <t>(y + r M )dy + ..., 

i.e., the approximation of the scaling function coefficient s° up to order M is 
made by sampling f(x) at the position t m + k. 

Note that all of the above calculations could have been carried out for 
the first derivative of f(x) giving an approximation to the scaling function 
coefficients, of f'(x): 

= fij + fy + /( M+1 )( r + k) J y M (j>{y + r)dy + .... 

It was assumed above that there exist one tm such that 

J <f>(x -f Tju )x m dx = 0, 

for m = 1, ..., M — 1. For m = 1 this tm is easy to find: 

J <f>(x + r M )xdx = I <f>{x)(x - T M )dx 

= J x<f>(x)dx — tm J <j>(x)dx. 

But / phi(x)dx = 1, therefore, 


tm 



That is, t m is simply the first moment of <t>{x). To find tm for m > 1 
the calculations are simple but a bit longer and require the result from the 
following theorem to show that there is one tm which is the same for all 
m = 1, ..., M — 1. 

If f 4>{x)dx = 1 and there exists r such that / <j>(x + T)x m dx = 0 for 
m = 1, ..., M — 1 then / <f>(x)x m dx = (/ <f>(x)xdx) m for m = 1, ..., M — 1. 
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Proof: Start with 

J <f>(x + r)x m dx = 0, 

and let y = x + r to get, 

J 4>{y)(y - T ) m = o. 

Using the binomial theorem this becomes, 

j <f>(y) E ^ ™ ) y r (- T ) m ~ rd y - °- 

Let the moments of <j>{x) be denoted by Mi = f cj)(x)x l dx to get 

E ( 7 ) M ' = 0 

A simple calculation yields r = M\. Using this value of r and summing only 
up to m — 1 the previous expression becomes, 

E ( 7) = 0. 

Or, 

= -E ( 7 ) 

From the hypotheses it is known that M 0 = f <f>{x)dx = 1. Therefore, M p = 
M 1 for p = 0,1, and with this knowledge it is easy to show that M p — 
for p = 2: 

= -E ( 7 ) 

which holds for m = 1,2. Combine the powers of M\ to get, 

= -wrE ( 7 ) 

But, this is nothing more than, 

M m = — Mj m [(l - l) m - 1], 
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or simply, 


Mm = Mr, 

where m = 1,2. The proof is complete, since higher powers of m can be 
found by induction. 

0.14 Appendix 2 

In this appendix the moments of <f>(x) will be calculated in closed form. Begin 
with the definition of the scaling function, 

<f>( x ) = - *)• 

k 

Next, calculating the m-th moment of <f>(x) yields, 

J <f>(x)x m — J <j>{ 2x — k)x m dx . 

Change the variable of integration such that y = 2x — k to get, 

/<H *)*" = E**/ *(y)(l/2)"(y + k) m l/2dy, 

k 

= (i/2r +i e ^ / m(s + •=)"<>»• 

k J 

Now, recall the binomial theorem to get, 

J <f>(x)x m = (l/2) m+1 ^2 h k J 4>{y) jr ^ ™ ^ y l k m ~ l dy 

Rewrite the moments of <f>(x) as Mi = / x t <f>(x)dx to get, 

M m = (l/2) m+1 £( 7) 

1=0 \ 1 / k 


Now let yi — Yd; hkk 1 to get 


m—l 


M t . 
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